Multi i-TED coincidence resolution¶

First¶

  • The current setup is considered the "final" multi i-TED system
  • Characteristics of the analysis:
    • Normalized
    • Background was not subtracted
    • Calibrated with 88c, 100ns
    • Thresholds: 888 and 88c
    • Integration windows: 100, 150, 200, 250 ns
    • Cs137 and Na22
  • Detectors
    • iTED A: Moving
    • iTED B,C,D: Closest
In [1]:
import panel
import ipywidgets

panel.extension('ipywidgets')
In [2]:
pkg_ver = lambda pkg: "{:<20}{:}".format(pkg.__name__,pkg.__version__)

# ROOT
import uproot
print(pkg_ver(uproot))
import ROOT

# Machine Learning
import sklearn
print(pkg_ver(sklearn))
import torch
print(pkg_ver(torch))

# Data science
import scipy
print(pkg_ver(scipy))
import numpy
print(pkg_ver(numpy))
import pandas
print(pkg_ver(pandas))
import awkward
print(pkg_ver(awkward))

# Visualizations
import matplotlib
print(pkg_ver(matplotlib))
import matplotlib.pyplot as plt
import hvplot
print(pkg_ver(hvplot))
import hvplot.pandas
pandas.options.plotting.backend = 'holoviews'

import tqdm
print(pkg_ver(tqdm))
from tqdm.auto import tqdm
tqdm.pandas()

import glob
uproot              5.0.7
Welcome to JupyROOT 6.28/02
sklearn             1.2.2
torch               2.0.0
scipy               1.10.1
numpy               1.23.5
pandas              1.5.3
awkward             2.2.0
matplotlib          3.7.1
hvplot              None
tqdm                4.62.3
In [3]:
%jsroot
In [4]:
"""
file = uproot.open(Data["Na22-40"](100))

c = file['COINCIDENCES;1'].arrays(library="ak")

# Gets ttree as dataframe
df = awkward.to_dataframe(c)
# Only multiplicity 2
df = df.query("detector_multiplicity == 2 & deposited_energy != 0")
# Only if it includes the scatterer of A
df = df.loc[df[df.index.get_level_values('subentry').isin([0])].index.get_level_values('entry')]
# Only if it only includes crystals of A
df = df[df.index.get_level_values('subentry').isin([0,1,2,3,4])]
# Only if 2 entries according to the previous filters
df = df.loc[df.index.get_level_values('entry')[df.index.get_level_values('entry').duplicated()]]
# Apply energy calibration to individual entries
df["deposited_energy"] = df.progress_apply(lambda x: numpy.polyval(calibrations_df.query(f"crystal == 'A{x.name[1]}' & cw == '100'").iloc[0][[0,1,2]][::-1],x.deposited_energy), axis=1)
# Recalculate sum of energies
new_sum = df.groupby(by="entry", axis=0).deposited_energy.sum()
df["total_deposited_energy"] = df.progress_apply(lambda x: new_sum.loc[x.name[0]], axis=1)
# Show dataframe
df

hist = ROOT.TH1F("","Coincidences;Energy [keV];Counts",2000,0,2000)
for i in df.iterrows():
    hist.Fill(i[1].total_deposited_energy)
    
TH1D_draw(hist).Draw()
"""
Out[4]:
'\nfile = uproot.open(Data["Na22-40"](100))\n\nc = file[\'COINCIDENCES;1\'].arrays(library="ak")\n\n# Gets ttree as dataframe\ndf = awkward.to_dataframe(c)\n# Only multiplicity 2\ndf = df.query("detector_multiplicity == 2 & deposited_energy != 0")\n# Only if it includes the scatterer of A\ndf = df.loc[df[df.index.get_level_values(\'subentry\').isin([0])].index.get_level_values(\'entry\')]\n# Only if it only includes crystals of A\ndf = df[df.index.get_level_values(\'subentry\').isin([0,1,2,3,4])]\n# Only if 2 entries according to the previous filters\ndf = df.loc[df.index.get_level_values(\'entry\')[df.index.get_level_values(\'entry\').duplicated()]]\n# Apply energy calibration to individual entries\ndf["deposited_energy"] = df.progress_apply(lambda x: numpy.polyval(calibrations_df.query(f"crystal == \'A{x.name[1]}\' & cw == \'100\'").iloc[0][[0,1,2]][::-1],x.deposited_energy), axis=1)\n# Recalculate sum of energies\nnew_sum = df.groupby(by="entry", axis=0).deposited_energy.sum()\ndf["total_deposited_energy"] = df.progress_apply(lambda x: new_sum.loc[x.name[0]], axis=1)\n# Show dataframe\ndf\n\nhist = ROOT.TH1F("","Coincidences;Energy [keV];Counts",2000,0,2000)\nfor i in df.iterrows():\n    hist.Fill(i[1].total_deposited_energy)\n    \nTH1D_draw(hist).Draw()\n'
In [5]:
class spectrum:
    
    def __init__(self, File_, Configuration_, Window_, Time_):  

        self.__FileName = File_
        self.__Configuration = Configuration_
        self.__Window = Window_
        self.__Time = Time_
        self.__File = uproot.open(File_)
        self.__DF = self.DF()
        self.__TH1D = self.TH1D()
        
    def __call__(self, ch):
        return numpy.polyval(self.__Calibration[::-1],ch)

    def File(self):
        return self.__File
    
    def FileName(self):
        return self.__FileName
    
    def DF(self):
        # Gets ak dataframe
        df_ak = self.__File['COINCIDENCES;1'].arrays(library="ak")
        #Transform from ak to pandas
        df = awkward.to_dataframe(df_ak)
        # Only multiplicity 2
        df = df.query("detector_multiplicity == 2 & deposited_energy != 0")
        # Only if it includes the scatterer of A
        df = df.loc[df[df.index.get_level_values('subentry').isin([0])].index.get_level_values('entry')]
        # Only if it only includes crystals of A
        df = df[df.index.get_level_values('subentry').isin([0,1,2,3,4])]
        # Only if 2 entries according to the previous filters
        df = df.loc[df.index.get_level_values('entry')[df.index.get_level_values('entry').duplicated()]]
        # Apply energy calibration to individual entries
        df["deposited_energy"] = df.progress_apply(lambda x: numpy.polyval(
            calibrations_df.query(f"crystal == 'A{x.name[1]}' & cw == '100'").iloc[0][[0,1,2]][::-1],x.deposited_energy
        ), axis=1)
        # Recalculate sum of energies
        new_sum = df.groupby(by="entry", axis=0).deposited_energy.sum()
        df["total_deposited_energy"] = df.progress_apply(lambda x: new_sum.loc[x.name[0]], axis=1)
                        
        return df
    
    def TH1D(self):
        # Create histogram
        hist_s = ROOT.TH1F("","Coincidences (Sum);Energy [keV];Counts",2000,0,2000)
        hist_1 = ROOT.TH1F("","Coincidences (Absorber 1);Energy [keV];Counts",2000,0,2000)
        hist_2 = ROOT.TH1F("","Coincidences (Absorber 2);Energy [keV];Counts",2000,0,2000)
        hist_3 = ROOT.TH1F("","Coincidences (Absorber 3);Energy [keV];Counts",2000,0,2000)
        hist_4 = ROOT.TH1F("","Coincidences (Absorber 4);Energy [keV];Counts",2000,0,2000)
        # Add entries
        for i in self.__DF.iterrows():
            match i[0][1]:
                case 1:
                     hist_1.Fill(i[1].total_deposited_energy)
                case 2:
                     hist_2.Fill(i[1].total_deposited_energy)
                case 3:
                     hist_3.Fill(i[1].total_deposited_energy)
                case 4:
                     hist_4.Fill(i[1].total_deposited_energy)
                case _:
                    continue
            hist_s.Fill(i[1].total_deposited_energy)
        # Scale
        hist_1.Scale(1/self.__Time)
        hist_2.Scale(1/self.__Time)
        hist_3.Scale(1/self.__Time)
        hist_4.Scale(1/self.__Time)
        hist_s.Scale(1/self.__Time)
        
        return [hist_s, hist_1, hist_2, hist_3, hist_4]
    
    def get_TH1D(self, i=0):        
        return self.__TH1D[i]
    
    def Configuration(self):
        return self.__Configuration
    
    def Window(self):
        return self.__Window
    
    def Rate(self):
        return self.get_TH1D().Integral()
In [6]:
calp = glob.glob('/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/Calibration/**/*.CALp', recursive=True)

calibrations = []

for i in calp:
    line = pandas.read_csv(
        i,
        sep = "      :  ",
        skiprows=[0,4,5,6],
        header=None,
        engine="python"
    ).drop([0], axis=1).T
    
    line["crystal"] = i.split("/")[-1][:2]
    line["cw"] = i.split("/")[-1][5:8]
    
    calibrations.append(line)
    
calibrations_df = pandas.concat(calibrations)
In [7]:
Data = {
    #"Na22-0":    lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_10_18_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x0.0_y0.0_CW{win}.root",
    #"Na22-10":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_25_35_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x10.0_y0.0_CW{win}.root",
    #"Na22-20":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_40_53_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x20.0_y0.0_CW{win}.root",
    #"Na22-30":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.18_56_16_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x30.0_y0.0_CW{win}.root",
    #"Na22-40":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_10/Na22_center_iTEDA_moving_D.2023_05_10_T.19_11_36_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x40.0_y0.0_CW{win}.root",
    #"Cs137-0":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.09_51_47_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x0.0_y0.0_CW{win}.root",
    #"Cs137-1":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_17/Cs137_mid_iTEDA_far_D.2023_05_17_T.20_57_08_C.itedABCD_lab_custom_2023.02.22_4.0v_887_1200s_serie2_2_CW{win}.root",
    "Cs137-0":   lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_18/Cs137_center_iTEDA_far_D.2023_05_18_T.18_31_25_C.itedABCD_lab_custom_2023.02.22_4.0v_887_1200s_CW{win}.root",
    "Cs137-10":  lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_07_01_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x10.0_y0.0_CW{win}.root",
    "Cs137-20":  lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_22_17_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x20.0_y0.0_CW{win}.root",
    "Cs137-30":  lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_37_33_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x30.0_y0.0_CW{win}.root",
    "Cs137-40":  lambda win: f"/run/media/bgameiro/d043b5e4-57a4-457e-8839-cb1adc9c72bc/Data/data_2023_05_11/Cs137_center_iTEDA_moving_D.2023_05_11_T.10_52_49_C.itedABCD_lab_custom_2023.02.22_4.0v_887_900s_x40.0_y0.0_CW{win}.root",
}
In [8]:
def get_resolution(spectr,en):
    
    TH1D = spectr.get_TH1D()
    
    TH1D.GetXaxis().SetRange(TH1D.FindBin(en-250),TH1D.FindBin(en+200))
    
    MaxBin   = TH1D.FindBin(TH1D.GetMaximumBin())
    
    ADC_Low  = MaxBin-150
    ADC_High = MaxBin+150

    gaussFit = ROOT.TF1("gaussFit", "gaus(0)+pol3(3)", ADC_Low, ADC_High)
    gaussFit.SetParameters(TH1D.GetMaximum(),MaxBin,10,0.375,-0.0005,0,0)
    TH1D.Fit(gaussFit,"QR")
    
    sigma = abs(gaussFit.GetParameter(2))
    centroid = gaussFit.GetParameter(1)
            
    return sigma*numpy.sqrt(2*numpy.log(2))*2/centroid*100, centroid, abs(centroid-en)/en*100
In [43]:
def TH1D_draw(spectr,en):
    
    canvas = ROOT.TCanvas()
    canvas.cd()
    
    TH1D = spectr.get_TH1D()
    
    #TH1D.SetTitle(repr(cell))
    TH1D.SetStats(False)
    
    latex = ROOT.TLatex()
    latex.SetNDC()
    latex.SetTextSize(0.03)
    
    TH1D.Draw("PE PLC")
    
    l1,l2,l3 = get_resolution(spectr,en)
        
    latex.DrawText(0.7, 0.8, "R: {:.2f}%".format(l1))
            
    latex.DrawText(0.7, 0.7, "E: {:.0f}keV".format(l2))
    
    TH1D1 = spectr.get_TH1D(1)
    TH1D1.Draw("SAME PLC")
    TH1D2 = spectr.get_TH1D(2)
    TH1D2.Draw("SAME PLC")
    TH1D3 = spectr.get_TH1D(3)
    TH1D3.Draw("SAME PLC")
    TH1D4 = spectr.get_TH1D(4)
    TH1D4.Draw("SAME PLC")
    
    legend = ROOT.TLegend(0.1,0.7,0.9,0.9)
    legend.AddEntry(TH1D,"Sum of absorbers","l")
    legend.AddEntry(TH1D1,"Absorber 1","l")
    legend.AddEntry(TH1D2,"Absorber 2","l")
    legend.AddEntry(TH1D3,"Absorber 3","l")
    legend.AddEntry(TH1D4,"Absorber 4","l")
    legend.Draw("SAME")
    
    canvas.BuildLegend()
            
    return canvas
In [10]:
entries = []

for i in tqdm(Data):
    source, conf = i.split("-")
    for CW in tqdm([100,150,200,250]):
        
        spectr = spectrum(
            Data[i](CW),
            int(conf), 
            CW,
            900 if "900s" in Data[i](CW) else 1200
        )
        
        en = 662 if source == "Cs137" else 511
        
        res, cent, fit = get_resolution(spectr,en)
                                    
        entries.append(
            pandas.DataFrame({
                "resolution": res,
                "fit": fit,
                "cps": spectr.Rate(),
                "cw": CW,
                "configuration": int(conf),
                "obj": spectr,
            }, index=[0])
        )
                    
entries_df = pandas.concat(entries, ignore_index=True)
  0%|          | 0/5 [00:00<?, ?it/s]
  0%|          | 0/4 [00:00<?, ?it/s]
  0%|          | 0/68550 [00:00<?, ?it/s]
  0%|          | 0/68550 [00:00<?, ?it/s]
  0%|          | 0/61970 [00:00<?, ?it/s]
  0%|          | 0/61970 [00:00<?, ?it/s]
  0%|          | 0/58044 [00:00<?, ?it/s]
  0%|          | 0/58044 [00:00<?, ?it/s]
  0%|          | 0/55510 [00:00<?, ?it/s]
  0%|          | 0/55510 [00:00<?, ?it/s]
  0%|          | 0/4 [00:00<?, ?it/s]
  0%|          | 0/52188 [00:00<?, ?it/s]
  0%|          | 0/52188 [00:00<?, ?it/s]
  0%|          | 0/46988 [00:00<?, ?it/s]
  0%|          | 0/46988 [00:00<?, ?it/s]
  0%|          | 0/43842 [00:00<?, ?it/s]
  0%|          | 0/43842 [00:00<?, ?it/s]
  0%|          | 0/42110 [00:00<?, ?it/s]
  0%|          | 0/42110 [00:00<?, ?it/s]
  0%|          | 0/4 [00:00<?, ?it/s]
  0%|          | 0/77886 [00:00<?, ?it/s]
  0%|          | 0/77886 [00:00<?, ?it/s]
  0%|          | 0/69772 [00:00<?, ?it/s]
  0%|          | 0/69772 [00:00<?, ?it/s]
  0%|          | 0/64626 [00:00<?, ?it/s]
  0%|          | 0/64626 [00:00<?, ?it/s]
  0%|          | 0/61522 [00:00<?, ?it/s]
  0%|          | 0/61522 [00:00<?, ?it/s]
  0%|          | 0/4 [00:00<?, ?it/s]
  0%|          | 0/97626 [00:00<?, ?it/s]
  0%|          | 0/97626 [00:00<?, ?it/s]
  0%|          | 0/87606 [00:00<?, ?it/s]
  0%|          | 0/87606 [00:00<?, ?it/s]
  0%|          | 0/81152 [00:00<?, ?it/s]
  0%|          | 0/81152 [00:00<?, ?it/s]
  0%|          | 0/77240 [00:00<?, ?it/s]
  0%|          | 0/77240 [00:00<?, ?it/s]
  0%|          | 0/4 [00:00<?, ?it/s]
  0%|          | 0/121068 [00:00<?, ?it/s]
  0%|          | 0/121068 [00:00<?, ?it/s]
  0%|          | 0/108498 [00:00<?, ?it/s]
  0%|          | 0/108498 [00:00<?, ?it/s]
  0%|          | 0/100654 [00:00<?, ?it/s]
  0%|          | 0/100654 [00:00<?, ?it/s]
  0%|          | 0/95422 [00:00<?, ?it/s]
  0%|          | 0/95422 [00:00<?, ?it/s]
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
In [11]:
entries_df["resolution"] = entries_df.obj.map(lambda x: get_resolution(x,662)[0])
entries_df["fit"] = entries_df.obj.map(lambda x: get_resolution(x,662)[2])
entries_df
Out[11]:
resolution fit cps cw configuration obj
0 8.940606 2.044403 19.550833 100 0 <__main__.spectrum object at 0x7f434d1fb820>
1 8.883976 1.959143 17.604167 150 0 <__main__.spectrum object at 0x7f434a2d3a60>
2 8.816424 1.898117 16.409167 200 0 <__main__.spectrum object at 0x7f434a275a20>
3 9.137346 1.949307 15.675833 250 0 <__main__.spectrum object at 0x7f434878a950>
4 8.636753 3.762054 19.307778 100 10 <__main__.spectrum object at 0x7f434a3eb520>
5 8.574016 3.721391 17.285556 150 10 <__main__.spectrum object at 0x7f43482b1c90>
6 8.491103 3.632212 16.070000 200 10 <__main__.spectrum object at 0x7f42d1f7ae30>
7 8.463257 3.594018 15.391111 250 10 <__main__.spectrum object at 0x7f434a3eb640>
8 8.926223 4.546599 29.617778 100 20 <__main__.spectrum object at 0x7f42d13ed5d0>
9 8.420145 4.330549 26.447778 150 20 <__main__.spectrum object at 0x7f433bfd8250>
10 8.378418 4.257623 24.390000 200 20 <__main__.spectrum object at 0x7f4349818d90>
11 8.314041 4.196396 23.146667 250 20 <__main__.spectrum object at 0x7f42d08597b0>
12 8.873971 5.154211 37.441111 100 30 <__main__.spectrum object at 0x7f42c9ee2b30>
13 8.794574 5.072420 33.440000 150 30 <__main__.spectrum object at 0x7f42d1972a70>
14 8.721396 4.970308 30.840000 200 30 <__main__.spectrum object at 0x7f42cf450820>
15 8.715796 4.923926 29.320000 250 30 <__main__.spectrum object at 0x7f42ce866350>
16 8.833005 4.799102 47.157778 100 40 <__main__.spectrum object at 0x7f42d0c39270>
17 8.804763 4.706591 42.051111 150 40 <__main__.spectrum object at 0x7f42d02ed4e0>
18 8.805299 4.643221 38.850000 200 40 <__main__.spectrum object at 0x7f434a17fc70>
19 8.819354 4.616517 36.714444 250 40 <__main__.spectrum object at 0x7f42d18722f0>
In [24]:
entries_df.hvplot.scatter("configuration","resolution",by="cw")
Out[24]:
In [13]:
entries_df.hvplot.scatter("configuration","fit",by="cw")
Out[13]:
In [14]:
entries_df.hvplot.scatter("configuration","cps",by="cw")
Out[14]:
In [44]:
l=entries_df.loc[0].obj
get_resolution(l,662)
TH1D_draw(l,662).Draw()
In [46]:
l=entries_df.loc[16].obj
get_resolution(l,662)
TH1D_draw(l,662).Draw()